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ABSTRACT 

We performed a series of three-dimensional numerical simulations of supersonic homogeneous Euler turbu- 
lence with adaptive mesh refinement (AMR) and effective grid resolution up to 1024 3 zones. Our experiments 
describe nonmagnetized driven supersonic turbulent flows with an isothermal equation of state. Mesh refine- 
ment on shocks and shear is implemented to cover dynamically important structures with the highest resolution 
subgrids and calibrated to match the turbulence statistics obtained from the equivalent uniform grid simulations. 

We found that at a level of resolution slightly below 512 3 , when a sufficient integral/dissipation scale separa- 
tion is first achieved, the fraction of the box volume covered by the AMR subgrids first becomes smaller than 
unity. At the higher AMR levels subgrids start covering smaller and smaller fraction of the whole volume that 
scale with the Reynolds number as Re' 1 ! 4 . We demonstrate the consistency of this scaling with a hypothesis 
that the most dynamically important structures in intermittent supersonic turbulence are strong shocks with a 
fractal dimension of two. We show that turbulence statistics derived from AMR simulations and simulations 
performed on uniform grids agree surprisingly well, even though only a fraction of the volume is covered 
by AMR subgrids. Based on these results, we discuss the signature of dissipative structures in the statistical 
properties of supersonic turbulence and their role in overall flow dynamics. 

Subject headings: ISM: structure — ISM: clouds — hydrodynamics — turbulence — methods: numerical 



1. INTRODUCTION 

Turbulence in molecular clouds is characterized by very 
high integral scale Reynolds numbers, Re = £()Uq/v > 10 8 , 
where £q « 10 pc is the typical scale on which the turbu- 
lence is driven, uq » 2 km s -1 is the typical velocity asso- 
ciated with that scale, and u is the kinematic viscosity of 
the molecular gas (see Elmegreen & Scalo 2004 for a re- 
view). Nonlinear interactions that dominate the dynamics of 
such multi-scale flows critically depend on adequate resolu- 
tion and suggest to exploit spatial and temporal adaptivity 
in numerical simulations. While adaptive mesh refinement 
has previously been applied to simulate the evolution of indi- 
vidual sing ular structur e s in incompressi ble Euler turbulence 
JPumir & Siggial Q990; iGrauer. Marliani. & Germ aschewski 
1 19981) . no attempts have been made so far to address the 
isotropic turbulence case with AMR. 

In this letter we apply high order adaptive methods to 
simulate supersonic turbulence in star forming molecular 
clouds with very high resolution. Such simulations are es- 
sential for understanding the nature of turbulent fragmenta- 
tion that may ultimately c ontrol the stellar initial mass func- 
tion JPadoan & Nordlundl2002l) . We exploit the fact that tur- 
bulent flows are not completely chaotic. Order is always 
present on both large and small scales due to the intermit- 
tent nature of turbulence. For instance, large under-dense 
voids and sharp density peaks are known to be characteris- 
tic of supersonic turbulent flows with a "soft" equation of 
state JPassot & Vazauez-Semadenilll998l) . We therefore ex- 
pect AMR simulations of isothermal turbulence to be bene- 
ficial in terms of computational resources. 1 We show that in 
fact AMR technique can be profitably applied to turbulence 
simulations in contrast with the prevailing belief that in such 

1 Application of AMR techniques to multiphase interstellar turbulence 
simulations, where the effective adiabatic index determined by the balance 
between heating and cooling at high densities falls below unity, appears to be 
even more promising. 



studies the adaptive mesh does not help as the fine structures 
emerge through the entire computational domain. 

2. TURBULENCE, INTERMITTENCY, AND ADAPTIVE MESHES 

It follows from the iKolmo goro vl d!941h phenomenology for 
turbulent cascade (K41) that the inertial range spans an in- 
terval of scales Af = lo/v ~ [V/(^o M o)] ~ Re 3 ! 4 , where 
V ~ {v 3 je) 1 ! 4 is the Kolmogorov dissipation scale, and the en- 
ergy cascade rate e ~ Uq/£q. Hence, the number of zones per 
integral scale £^ required to model a three-dimensional tur- 
bulent flow using a finite-difference method is A/" 3 oc Re 9 / 4 . 
The storage resource for a numerical experiment of this sort 
would also grow as S oc Re 9 ! 4 , while the number of CPU 
hours needed to simulate high-Reynolds-number flows on a 
uniform grid for a given number of large eddy turnover times 
would scale as T oc Re 3 (e.g. Frisch 1995). 

However, the above calculation implicitly assumes that the 
inertial range flow is completely chaotic and does not con- 
tain any coherent structures, which is often not the case in 
experimental \\\g\\-Re flows. In fact, both laboratory ex- 
periments and numerical simulations indicate the presence 
of some order in the flows on both small I ~ rj and large 
I ~ £q scales. Since turbulence is intermittent, the K41 the- 
ory gives only an approximation to its statistical properties 
and requires corrections to reproduce experimental measure- 
ments of high-order statistics. One way to add a form of 
intermitten cy to the K41 mod el is known as the /3-model 
(Frisch, Sulem, & Nelkin 1978). It assumes that the fraction 
of the space occupied by the 'active' eddies of size I scales 
as (£/£q) 3 ~ d , where D is interpreted as the fractal dimension 
of small-scale dissipative structure. In the framework of the 
/3-model the dissipation scale is a function of D, therefore 
JV D = £o/v(P) ~ Re 3 ^ 1+D) . One can easily recover the fa- 
miliar K41 result from this formula by setting D = 3 which 
corresponds to the zero-intermittency limit. Thus, the scaling 
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FIG. 1 . — Gas density from AMR simulation of supersonic Mach 6 turbulence with effective resolution of 1024 3 at time t = 6 tj y „. Left panel: a projection of 
the density field through the computational box (log scale, dense regions are dark). Right panel: a thin slice from the same density distribution with contours 
showing the patches of the base grid not covered by the AMR subgrids. The total volume covering factor of the first level AMR subgrids is approximately 65%. 



for the storage depends on the fractal dimension of the small- 
scale dissipative structure, Sq oc Re 3D ^ 1+D \ 

There are good reasons to assume D w 1 in the incom- 
pressible case, where vortex filaments are the most si ngular 
dissipative structures (She & Leveaue 1994; Dubrulle 1994), 
and D w 2 for supersonic compressible turbulence, where the 
dissipative structures are shocks (Boldvrev 2002). For these 
fractal dimensions, Si /S ~ Re~ 3 / 4 and S2/S ~ Re~ l l A , so the 
advantage of adaptive methods for simulations of high-/?<? tur- 
bulent flows is quite clear. 

There is a simple direct analogy between the process of 
building an AMR hierarchy to resolve the turbulent structures 
of the fractal dimension D and the box-counting method for 
determination of that same fractal dimension. The fractal di- 
mension is given by the relation n(S) oc S~ D , where n{5) is the 
number of boxes of linear size 5 needed to cover the set of sin- 
gular structures. One can think of the box size S as a function 
of the level number / in the AMR hierarchy, Si = Sof~ l , where 
So is the zone size on the base grid and / is the refinement fac- 
tor (usually / = 2 or 4). Then the number of zones needed to 
cover the structures on level I would grow as «/ oc Sj D oc f ID , 
and the volume covering factor of level I subgrids would scale 
with the level number as F[ oc f l< - D ~ 3 \ Assuming D = 2, one 
gets F{ oc f~ l . Having chosen an appropriate resolution for the 
base grid and / = 4, one can predict the covering factors for 
/ = 1, 2, and 3-subgrids to be 0.25, 0.06, and 0.015, respec- 
tively. For incompressible turbulence, the same factors are 
valid iff = 2. 

At what base grid resolution AMR first becomes benefi- 
cial for turbulence simulations? As with the box-counting 
method, the expected scaling can be attained only with a 
high enough resolution. It is necessary to provide inte- 
gral/dissipation scale separation for the inertial range to be- 
come extended in the wavenumber space. In simulations of 
compressible Euler turbulence with PPM this can be achieved 
at resolutions somewhere in between 256 3 and 512 3 (schemes 
with higher numerical diffusivity than PPM would require 
many more zones), while for compressible Navier-Stokes tur- 



bulence one needs about 4 3 times higher resolution since a 
larger range of scales is needed to account for the physical 
dissipation (Svtineetal. 2000). Therefore, in PPM simula- 
tions of Euler turbulence the advantage of AMR comes into 
play earlier than in those of Navier-Stokes turbulence. On the 
other hand, due to the difference in scaling behavior, AMR 
is expected to save more resources in numerical experiments 
with incompressible flows at very high Reynolds numbers. In 
the following sections we show that in practice for Mach 6 
isothermal turbulence the advantage of AMR is indeed first 
noticeable at a resolution of about 512 3 , and the scaling rela- 
tions derived above are consistent wi th our numerical exp eri- 
ments which involve the Enzo code (O'Shea et al. 2004J and 
references therein). 

3. NUMERICAL APPROACH 

The Enzo code uses a dire ct Eulerian formulation of 
the Piecewise Parabolic Method ( Colella & Woodward 1984, 
PPM) to solve the equations o f gas dynamics on a hie rarchy 
of structured adaptive meshes (Berger & Colella 1989). 

In order to mimic the conditions in molecular clouds, we 
adopt a quasi-isothermal equation of state with the ratio of 
specific heats 7 = 1.001. To test our PPM implementation in 
this highly compressible low-7 regime and at high Mach num- 
bers, we ran a number of tests with and without AMR, includ- 
ing simple shock tube tests described in (Balsara 1994). We 
found very good agreement with the exact solutions and with 
results obtained with a Riemann solver for purely isothermal 
gas dynamics. The temperature variations in tests with Mach 
6 flows, due to the small deviation of 7 from unity, were typi- 
cally below 0.1%. 

The AMR subgrids are placed where needed to resolve 
shocks with large pressure jumps and regions of strong shear. 
We identify shocks using the PPM shock detection algorithm. 
In addition, we also use a norm of the velocity gradient ma- 
trix || diiij || to account for shear. This matrix norm is simi- 
lar to the Frobenius norm, except that it does not include the 
contribution from diagonal elements. To get a better AMR 
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FIG. 2. — Volumetric rendering of the gas density in a 100-zone-thick slice 
through a subsample of the computational domain. V- and U-shaped shock- 
lets or "Mach cones" are the most common structures in supersonic isother- 
mal hydrodynamic turbulence. As in a hierarchy of vortices in incompressible 
turbulence, large-scale Mach cones are broken into smaller and smaller ones 
down to the finest resolved scale. Mach angles of the largest cones visible in 
this snapshot correspond to relative fluid velocities of (2 — 3) c s . 




-3 



-1 



1 





I°910 P 

FIG. 3 . — Comparison for turbulence statistics based on 256 3 and 5 12 3 uni- 
form grid simulations and an AMR simulation with effective grid resolution 
of 1024 3 at time t = 6 tj y „. Probability density functions (PDFs) for the gas 
density. 

efficiency for simulations of turbulent flows (which involve a 
large number of subgrids) we use mesh refinement by a factor 
of 4. We obtained a very good agreement between our non- 
adaptive and AMR runs for refinement on shocks with pres- 
sure jumps Ap/p > 2. When a combination of shocks and 
shear controlled the refinement, a minimum pressure jump 
of 3 was sufficient, while the application of the second re- 
finement criterion accounted for some 20% more zones to be 
flagged. Note that these values are given only for orientation, 
while the actual lower bound for pressure jumps in refined 
shocks would depend on the AMR implementation. 

To maintain the turbulent kinetic energy in the computa- 
tional box at a given level, we use large-scale solenoidal force 
per unit mass with a fixed spatial pattern and a constant power 
in the range of wavenumbers k £ [1,2]. 




1.5 2 
log 10 k 

FIG. 4. — Same as in Fig. [3] but for the gas density power spectra. 
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FIG. 5. — Same as in Fig. |3] but for the compensated total velocity power 
spectra, scaled so that the curves match in the numerical dissipation range of 
PPM where they are invariant. 

4. TURBULENT STRUCTURES AT HIGH RESOLUTION 

We first performed a series of nonadaptive simulations of 
driven turbulence with an rms Mach number of 6 varying the 
uniform grid resolution from 128 3 up to 512 3 zones. We 
started the simulations with a uniform gas density and ran 
them for 6 dynamical times 2 to follow the development and 
saturation of the turbulent flow. We then restarted the 128 3 
run as an AMR simulation with the base grid of 128 3 and one 
level of refinement by a factor of 4, ran it from 4.8 to 6 td yn , 
and compared the results with our uniform grid run with the 
equivalent resolution of 512 3 zones. After a few iterations we 
found that refinement on shocks with pressure jumps above 
k 2 gives us velocity power spectra consistent with the spec- 
tra from nonadaptive run with the same effective resolution. 
The velocity statistics appear to be more sensitive to the re- 
finement criterion then the density PDF or the density power 
spectrum. Since the volume covering factor of the first level 
subgrids in this AMR simulation was about 90%, it remained 
unclear whether it is the high covering fraction that provides 
convergence or it is indeed the right refinement criterion. Be- 
fore switching to higher resolution base grids (it is wasteful 
to run AMR covering 90% of the volume) we restarted the 
same AMR run from t = 6 td yn allowing for two levels and got 
an estimate of the covering fraction for the second level about 

2 tdyn = L/(2M), where L is the box size, M is the rms Mach number, and 
the sound speed c s is unity. 
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34% at effective resolution of 2048 3 . This number should be 
considered as a lower limit since it will grow by a few percent 
over the following f,/ v „, while the flow gets fully resolved on 
the second level. 

We then repeated the same experiment with the 256 3 base 
grid and one level of refinement, with however a slightly dif- 
ferent refinement recipe which also included shear, see Sec- 
tion 3 for details. With effective resolution of 1024 3 this AMR 
run is the largest simulation to date of supersonic turbulence 
in molecular clouds. Fig. ^shows a snapshot of the gas den- 
sity field from this run. From the left panel, where we show a 
distribution of the column density, it may indeed look like the 
fine structures emerge through the entire computational do- 
main. However, the morphology of the projected density field 
is quite different from what one can see in a slice, highlight- 
ing a significant loss of information built into the projection 
procedure. As one can see from a thin slice through the box 
shown in the right panel, AMR does help in this case since 
turbulence is very intermittent. The first level subgrids do not 
cover the under-dense voids and also some smaller windows 
in the high and intermediate density gas where the strong 
shocks are absent. These regions do not contribute much to 
cascading the energy from large to small scales as their share 
in the turbulence statistics is the same whether they are refined 
or not. Strong shock interactions and associated nonlinear in- 
stabilities create a very sophisticated multiscale pattern in the 
dynamically active regions (Figs. HI3 that is morphologically 
similar to what is observed in molecular clouds. This pattern 
is missing in numerical simulations at lower Re. We identify 
Mach cones and U-shaped shocklets as self-similar elements 
of this pattern (Fig. [2j. These effects of nonlinear shear in- 
stabilities resolved in our simulations can explain the lack of 
large-s cale shock si gnatures in the observations of molecular 
gas bv lBruntl ll2"b03). 

While we do not have an equivalent nonadaptive simula- 
tion, the comparison with PDFs and power spectra from 256 3 
and 512 3 runs shown in Figs. [3]-0speaks for itself. The den- 
sity PDFs can be very precisely fitted by a log-normal distri- 
bution and the AMR data match those from the two nonadap- 
tive runs, although the superior resolution of the AMR run 
provides better sampling for the high end of the density dis- 
tribution. Interactions of strong counter-propagating shocks 
are responsible for intermittent oscillations in the high density 
wing of the PDF, which are well resolved in our AMR simu- 
lation. These interactions also cause transient strong rarefac- 
tions lagging behind in time. As a result, on time-scles short 
compared to the dynamical time, the density PDF slightly 
wanders around its average log-normal representation at the 
highest densities and displays large temporal oscillations in 
its low-density end, see Fig. [3] The same processes reveals 
itself in correlated variations of the three-dimensional power 
spectrum of density in the inertial range. The slope lies some- 
where between -0.8 and -0.9 in our nonmagnetized models, 
see Fig. |4] The spectrum gets shallower upon the collisions of 
strong shocks, when the PDF's high density wing rises above 



the average log-normal representation. 

The density power spectrum builds up quickly at high 
wavenumbers, after we switch the AMR machinery on, sim- 
ply because PPM starts resolving the shocks better. However, 
the relaxation of the velocity power proceeds slower since it 
takes about one dynamical time for the resolved local sharp 
density structures to get through nonlocal dynamical interac- 
tions involving multiple scales. Only then the inertial range 
really extends to smaller scales. When AMR is first activated, 
the velocity power is insufficient at k > 25 and scales approx- 
imately as k~ 2 in this range. It then steadily accumulates at 
those frequencies for about td yn and saturates exactly at the 
level predicted by our nonadaptive simulations, see Fig. [5] 
The slope of the velocity power spectrum for the snapshot 
shown is about -1 .85, i.e. somewh at steeper than a slope of 
-1 .74 predicted by Boldvrev (2002) for supersonic turbulence 
assuming D = 2 and assuming the third order structure func- 
tion exponent equal to unity, as in incompressible turbu lenc e. 3 
According to the formalism developed by Dubrulle ( 1994), a 
slope of -1.85 would imply D w 2.3. This is consistent with 
the dimensionality of dissipative structures that can be inde- 
pendently estimated using the volume covering fractions of 
AMR subgrids at different levels of resolution (90, 65, and 
34% at 512 3 , 1024 3 , and 2048 3 , respectively). Such estimate 
returns the same value of D w 2.3. Within the uncertainties, 
both independent estimates agree with the observationally de- 
termined fractal dimension of molecular clouds D = 2.3 ±0.3 
JElmegreen & Falgaronell99fl). 

5. CONCLUSIONS 

While details of AMR implementation may vary and may 
have to be further refined to reproduce higher order statistics, 
it is clear that adaptivity in both space and time is indispens- 
able for numerical experiments with homogeneous isotropic 
supersonic turbulence at very high Reynolds numbers with 
grid resolutions > 1024 3 . 

Our simulations reveal a pattern of small-scale structures 
that is completely missing at lower Re. These structures orig- 
inate in nonlinear instabilities inherent in isothermal super- 
sonic flows and may control the scaling properties of turbu- 
lence and the fractal dimension of the dynamically important 
structures. Since the presence of magnetic fields can modify 
the unstable modes, it cannot be taken for granted that mag- 
netized turbulence should have the same scaling properties. 

Numerical experiments at very high Reynolds numbers are 
crucial for studies of turbulence in molecular clouds. 

This work was partially supported by NRAC allocation 
MCA098020S and utilized computing resources provided by 
the San Diego Supercomputer Center. 

3 Averaging over independent realizations of the turbulent flow is needed 
to obtain reliable estimates for the scaling exponents. This lies outside the 
scope of this letter which is primarily focused on the applicability of adaptive 
methods. 
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